function MertonDensity
clc;close all;

% Model parameters

T         = 5;
r     = 0.05;
S0    = 100;
sigma = 0.2;
muJ   = 0;
xiP   = 1;

% Define the COS method integration range

L = 8;
a = - L * sqrt(T);
b = + L * sqrt(T);

% Domain for the density f(x)

y = linspace(0.0001,4*S0,70);
dy = y(2)-y(1);
x=log(y);

% Effect of sigmaJ

sigmaJV = [0.0, 0.15, 0.3, 0.45];
argLegend = cell(4,1);
idx = 1;
pdf_M = zeros(length(x),length(sigmaJV));
cdf_M = pdf_M;
for i=1:length(sigmaJV)
    sigmaJTemp = sigmaJV(i);

    % Compute ChF for the Merton model

    cf = ChFForMertonModel(r,T,muJ,sigmaJTemp,sigma,xiP,S0);

    % Reference 

    f_XExact = COSDensity(cf,x,2^14,a,b);    
    pdf_M(:,i)=1./y.*f_XExact;
    trapz(y,pdf_M(:,i))
    cdf_M(:,i)=cumsum(pdf_M(:,i).*dy);
    argLegend{idx} = sprintf('sigma_J=%.2f',sigmaJTemp);
    idx = idx + 1;
end
MakeFigure(y, pdf_M,argLegend,'','S(T)','PDF')
xlim([min(y),max(y)])

MakeFigure(y, cdf_M,argLegend,'','S(T)','CDF')
xlim([min(y),max(y)])
 
function cf = ChFForMertonModel(r,tau,muJ,sigmaJ,sigma,xiP,S0)

% Term for E(exp(J)-1)

i = complex(0,1);
helpExp = exp(muJ + 0.5 * sigmaJ * sigmaJ) - 1.0;
  
% Characteristic function for the Merton model    

cf = @(u) exp(i*u*log(S0)).*exp(i * u .* (r - xiP * helpExp - 0.5 * sigma * sigma) *tau...
        - 0.5 * sigma * sigma * u.^2 * tau + xiP * tau .* ...
        (exp(i * u * muJ - 0.5 * sigmaJ * sigmaJ * u.^2)-1.0));

function f_X = COSDensity(cf,x,N,a,b)
i = complex(0,1); %assigning i=sqrt(-1)
k = 0:N-1; 
u = k * pi / (b - a);

% F_k coefficients

F_k    = 2 / (b - a) * real(cf(u) .* exp(-i * u * a));
F_k(1) = F_k(1) * 0.5; % Multiply the first term

% Final calculation

f_X = F_k * cos(u' * (x - a));

function MakeFigure(X1, YMatrix1, argLegend,titleIn,xLab,yLab)

%CREATEFIGURE(X1,YMATRIX1)
%  X1:  vector of x data
%  YMATRIX1:  matrix of y data

%  Auto-generated by MATLAB on 16-Jan-2012 15:26:40

% Create figure

figure1 = figure('InvertHardcopy','off',...
    'Colormap',[0.061875 0.061875 0.061875;0.06875 0.06875 0.06875;0.075625 0.075625 0.075625;0.0825 0.0825 0.0825;0.089375 0.089375 0.089375;0.09625 0.09625 0.09625;0.103125 0.103125 0.103125;0.11 0.11 0.11;0.146875 0.146875 0.146875;0.18375 0.18375 0.18375;0.220625 0.220625 0.220625;0.2575 0.2575 0.2575;0.294375 0.294375 0.294375;0.33125 0.33125 0.33125;0.368125 0.368125 0.368125;0.405 0.405 0.405;0.441875 0.441875 0.441875;0.47875 0.47875 0.47875;0.515625 0.515625 0.515625;0.5525 0.5525 0.5525;0.589375 0.589375 0.589375;0.62625 0.62625 0.62625;0.663125 0.663125 0.663125;0.7 0.7 0.7;0.711875 0.711875 0.711875;0.72375 0.72375 0.72375;0.735625 0.735625 0.735625;0.7475 0.7475 0.7475;0.759375 0.759375 0.759375;0.77125 0.77125 0.77125;0.783125 0.783125 0.783125;0.795 0.795 0.795;0.806875 0.806875 0.806875;0.81875 0.81875 0.81875;0.830625 0.830625 0.830625;0.8425 0.8425 0.8425;0.854375 0.854375 0.854375;0.86625 0.86625 0.86625;0.878125 0.878125 0.878125;0.89 0.89 0.89;0.853125 0.853125 0.853125;0.81625 0.81625 0.81625;0.779375 0.779375 0.779375;0.7425 0.7425 0.7425;0.705625 0.705625 0.705625;0.66875 0.66875 0.66875;0.631875 0.631875 0.631875;0.595 0.595 0.595;0.558125 0.558125 0.558125;0.52125 0.52125 0.52125;0.484375 0.484375 0.484375;0.4475 0.4475 0.4475;0.410625 0.410625 0.410625;0.37375 0.37375 0.37375;0.336875 0.336875 0.336875;0.3 0.3 0.3;0.28125 0.28125 0.28125;0.2625 0.2625 0.2625;0.24375 0.24375 0.24375;0.225 0.225 0.225;0.20625 0.20625 0.20625;0.1875 0.1875 0.1875;0.16875 0.16875 0.16875;0.15 0.15 0.15],...
    'Color',[1 1 1]);

% Create axes

%axes1 = axes('Parent',figure1,'Color',[1 1 1]);
axes1 = axes('Parent',figure1);
grid on

% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[45 160]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes1,[19 26]);
% Uncomment the following line to preserve the Z-limits of the axes
% zlim(axes1,[-1 1]);

box(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
% plot1 = plot(X1,YMatrix1,'Parent',axes1,'MarkerEdgeColor',[0 0 0],...
%     'LineWidth',1,...
%     'Color',[0 0 0]);

plot1 = plot(X1,YMatrix1,'Parent',axes1,...
    'LineWidth',1.5);
set(plot1(1),'Marker','diamond','DisplayName',argLegend{1});
set(plot1(2),'Marker','square','LineStyle','-.',...
    'DisplayName',argLegend{2});
set(plot1(3),'Marker','o','LineStyle','-.','DisplayName',argLegend{3});
set(plot1(4),'DisplayName',argLegend{4});

% Create xlabel

xlabel(xLab);

% Create ylabel

ylabel(yLab);

% Create title

title(titleIn);

% Create legend

legend1 = legend(axes1,'show');
set(legend1,'Color',[1 1 1]);
